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TABLES OF SPHERICAL BESSEL FUNCTIONS OF THE FIRST AND SECOND 
KINDS FOR COMPLEX ARGUMENT 


I. INTRODUCTION 


Spherical Bessel functions arise in wave propagation problems of 
spherical geometry. When loss mechanismg are included, the arguments of these 
functions are complex. Although researchers requiring these functions 
generally develop their own computer routines, a need nonetheless exists for 
tabulated values to verify the accuracy of such routines. Due to the limited 
avallability of tabulated values of spherical Bessel functions of complex 
argument, at least one group of researchers [1] had to resort to re-expressing 
these functions in terms of hypergeowetric functions in order to proceed with 
the analysis of their theoretical expressions. The author of this report 
required values of these functions while working ou a similar theoretical 
problem [2]. To the author's knowledge, only one other table of spherical 
Bessel functions of complex argument has been published [3]. This previous 
table is not completely adequate for those wishing to verify a computer 
routine. This is because to properly verify this type of routine, four 
extreme cases of argument are required: 1) a large real and a large imaginary 
part, 2) a large real and a small imaginary part, 3) a small ceal and a large 
imaginary part, and 4) a small real and a small imaginary part (all these 
statements are relative to unity). Case 3 is absent in these previous 
tables. Also, due to the relatively large number of variable parameters in 
selecting arguments to include’ in such a table, there are significant gaps 
between the arguments that are tabulated there. The present tubles give 
values of spherical Bessel functions of the first two kinds for all integral 
combinations of real and imaginary arguments ranging from 1 to 10 (in addition 
to listing several extreme cases). The method described in Ref. 3 for 
computing the functions was based on a continued fraction schete. ‘The 
approach taken here is via use of cecurrence relations. Finally, no reliable 
method of verifying the tabulated values given in Ref. 3 was desuribed. 
Several different such verification methods are herein presented. At least 
one of these methods is directly applicable in a straightforward way to any 


routine that computes spherical Bessel functions. Hence, via use of such a 
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verification method, researchers may verify routines that compute values out 
of the range of the available tables. 


Section II describes the methods used to compute the functions 
themselves. Section III presents the method of verifying the accuracy of the 
tabulated values. Appendix A gives tables of the spherical Bessel functions 
of the first and second kinds, jy) and ¥n(2)> respectively, with the real 
and imaginary parts of the complex variable z each varying from 1 to 10 and 
with the order varying from 0 to 20. Appendix B presents 5, (2) and y,(z) for 


some more extreme values of order and argument. 


These tables are not intended to be extensive since their primary purpose 
is to help establish validity of computer routines. Nonetheless, the tables 
are sufficiently comprehensive so that at least approximate interpolation is 
possible for real and/or imaginary arguments in the range 1 to 10 and for 


orders between 0 and 20. 


II. METHODS OF COMPUTING THE FUNCTIONS 


The recurrence relations for spherical Bessel functions [4] formed the 


basis for their computation 


£,(2) + £24,(2) = (2ntl)2 f(z) Q) 


where £ 62) = either 5,02) or Ynl2)> the spherical Bessel functions of the 
first and second kinds, respectively. However, this kind of calculation 
contains several potential sources of difficulty and has to be handled 


separately for the cases of the functions of the first and second kinds*. 


In the case of the cylindrical Bessel functions of the first kind of real 
argument, tt has been shown that the recurrence relation must actually be 


iterated backward [5]. The reason for this is the fact that as the order n of 


-— = ee eee 


*There are many subtleties involved in using recurrence relations to analyze 
functions. For a more rigorous analysis than the one given here, the reader 
should consult a standard reference on numerical computation, such as the one 
by Luke [4a]. 
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JAG) increases for fixed argument x, the functional value approaches zero 


extremely rapidly. This can readily be seen from the expression [6] 


thee 
100 + = ($3) (2) 


(for v*~ through real positive values). 


Since these functions J,,(x) decrease with increasing order v, the 
repeated use of the appropriate recurrence relation with increasing values of 
the order leads to increasingly less accurate values of the functions due to 
subtraction errors. This difficulty is usually circumvented by requiring 
J,(%) =a, (where m>> |x| and a is a smali number relative to 1) and by working 
Eq. (1) backward to Jo(z)- In this case the subtraction error effect is 
reversed and the accuracy of each function is improved with each decrease in 
order (until the limit of the number of digits available on the digital 
computer being used is reached). When this process is completed, all the 
values of J,(z) are off by a common factor, providing m is chosen sufficiently 
larger than the highest desired order. This common factor can be computed for 
the case of cylindrical Bessel functions via use of an addition theorem [7], 


and then the computed functions can be properly normalized. 


This procedure can also be applied, with minor modifications, to the 


computation of the spherical Bessel functions of the first and second kinds of 
complex arguments, 5,2) and Yn(2)> The procedure is not restricted to real 
arguments (or even to Bessel functions, for that matter). The only essential 
requirement is that the functional values decrease systematically with 


decreasing or with increasing order. If the functions decrease with 


increasing order, the recurrence relation must be iterated backward, and vice 
versa. 

Consider first the case of 5, (2) It is related to the cylindrical 
Bessel function J\(z) via the relation [8] 


5,2) * Viz Jat Yf?)> 


Equations (2) and (3) may readily be combined to give 


(3) 
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(for ne through real positive values). 


(4) 


It is clear from Eq. (4) that j,(z) approaches zero when the quantity 
{ez/(2n+l)] becomes less than. 1. This certainly will occur when 


n> Sh, (5) 
When |z| is large (greater than about 50 or so), Eq. (5) may be used to 
establish a criterion for the minimum value of n required to start the reverse 
iteration of Eq. (1). Such a criterion must also, of course, take into 
account the number of decimal digits used to represent the real imaginary 
parts of complex numbers in the available digital computer. For example, if 
the computer reliably represents seven decimal digits, m must be chosen such 
that the magnitude of j,.(z) is at least eight orders of magnitude smaller than 
the magnitude of j,(2) for the highest desired order n. This ensures that the 
“reverse subtraction error” effect described above will result in seven 
veliable significant digits (with the last digit rounded). For arguments 
|z|<50 it is safer to simply start the recurrence relation from a fixed large 
value of n.- In the program used to generate the tables in this report this 
large fixed value of n was chosen to be 250 (which is well beyond the eight 


orders of magnitude criterion mentioned above). 


Although the above analysis demonstrates the fact that 4, (2) approaches 
zero for fixed z as n?™, it does not follow that this occurs in a uniform 
fashion. In fact, aside from the normal trigonometric oscillation associated 
with spherical Bessel functions, it is possible that the values of jy (2) 
increase with increasing z [and vice versa with y,(z)], up to some finite 


value of n prior to exhibiting their asymptotic behavior. A digit of accuracy 


is lost for each order of magnitude of change in 5,62) or ¥n 2) which is 
contrary to their asymptotic behavior with respect to order. Hence, it is 
necessary to perform the iterations of the recurrence relations in double 
precision arithmetic to assure single precision accuracy of the final 
results. (It is still possible to lose single precision accuracy if this 
contrary behavior continues for more than seven orders of magnitude; hence, a 


test procedure is required. This is discussed further in Section III.) 
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Another problem can arise if the initializing value a is chosen too large 
{. or if the iteration process must be carried through too many orders. In such 
: cases, since the functions j,(z) will grow with decreasing n, it is possible 
oe for the un-normalized functions to exceed the largest numerical value that can 

a be represented in the digital computer being used for the computation. This 

; difficulty can be circumvented by not using the relation (1) directly, but 

‘ rather a modified form of it. This modification is based on the fact that the 
ratio of two successive spherical Bessel functions approaches zero at a much 

a . lower rate than the functions themselves. This may be seen by using Eq. (4) 


( to analyze the ratio of two successive spherical Bessel functions r,(2): 


. Sati (2) yor 1/2) Qn+l n+(1/2) 
S Fat?) = FG wt (xas3)(*S2 (6) 


ae By use of the theorem [7a} > il + =)" = e”, it is straightforward to 


show that r?2/2n, as n?@. seacese Th approaches zero at the rate 1/n, while 
{ it can be shown (via the same limit theorem) that = approaches zero at the 
rate 1/{n(2n)"]. 


Hence, the problem of exceeding the limit of the greatest number 
representable on the available digital computer can be virtually eliminated by 


using the recurrence relation for r,(z): 


+ roe (2nt+1)/z. (7) 


Equation (7) is obtained from Eq. (1) by dividing by f(z) and using the 
definition r, = £.4,(2)/£,(z)- 


Calculating ratios via Eq. (7) instead of directly computing the 
functions themselves has an additional advantage: Since all the desired 
functions in the current method are off by a common factor, the ratios 
obtained by the method outlined above are all automatically correct. Hence, 
to obtain the correct final values for the functions, it only remains to 
compute one of the desired functions accurately via an independent 


technique. Once this function Ls known, all other orders can be computed 
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using these ratios. For example, if 5,2) is known, then 4, (2) can be 
| computed from T(z) by its definition 


2 4,(2) = £5(2)4,(2)- (8) 
This process can clearly be repeated to obtain jo(2), j3(2), etc. 


As straightforward as this process is, a potential hazard exists if the 
“known” functional value is carelessly chosen. In order to understand this, 
consider first the analytical expressions for the first two spherical Bessel 
functions of the first kind [9]: 


sinz 
j,(2) = 
sinz cosz 
gE gt gr (9) 


If j,(z) should be chosen as the function used for the normalization 
process, and if the real part of z happens to be close to zero, the functional 
values returned by the standard complex sine function of FORTRAN can be rather 
inaccurate [resulting in inaccurate values for all orders of 5,(2)]- This 
problem can be circumvented by either writing a subroutine to calculate sinz 
more accurately (using its elementary infinite series representation) or by 
normalizing via 3, (2) instead. In generating the tables listed in this report 
5, (2) was used for all required normalizations with an accurate subroutine 


used for sinz in the case of small real argument. 


One final difficulty arises during this normalization procedure. This 
difficulty occurs when the real part of the argument z is large relative to 


20. To understand this, consider the following expression for sinz: 


ei(xtiy)_-i(xtly) 
sinz = sin(xtly) = Seg (i0) 


7 
a 


a 


{e 7 (cosxtisinx) - e” (cosx-isinx)]}. 


When the real part of z [represented by x in Eq. (10)] becomes large 
relative to 2m, the value of the trigonometric functions calculated by the 
standard FORTRAN routines can no longer be regarded to be accurate. To avoid 
this problem in generating the tables herein provided, the value of x was 
always taken to be the smallest possible quantity in excess of the largest 
possible integral multiple of 21 which does not exceed x. (It is assumed that 
when large values of x are used, these values are accurate to a sufficient 
number of digits to justify this procedure.) This process can be seen to be 


valid by examination of the theorems: 


sin(wt27m) = sinw 
(11) 


cos(wt2mm) = cosw , 


where wt21m = x, and m is chosen to be the largest integer which 
simultaneously satisfies this definition and yields the same mathematical sign 


for both w and x. 


Essentially all of the above comments made regarding 4,2) are also 
applicable to the spherical Neumann function Yn)» with the exception that in 
this case the recurrence relation [Eq- (1)] must be iterated forward. This 
once again is the consequence of the fact that to avoid subtraction errors the 
recurrence relation must be used in the direction that causes the functional 


values to increase. In analogy to Eq. (2), we have for ¥,62) [10]: 
me 
yq(2) = yixe Yat yf?) (12) 
and in analogy to Eq. (3) we have [6] 


Y (2) + -¥2 (sz) (13) 


(for ve through real positive values). 


As before, Eqs. (12) and (13) may be combined to yield 


n+ l/, 
y (z) -~ = Sacra eons Se os (2224) (14) 
n ez 
y+ Vy )z 
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(for n> through real positive values). 


Equation (14) clearly shows that y,(z) + @ as n + ™, requiring forward 
iteration of recurrence relation. As with 5,2), however, it is safer to use 
the recurrence relation of ratios in order to avoid the problem of exceeding 
the larger number representable by the available digital computer. 
Difficulties similar to those encountered in analyzing j,(2) arise during the 
normalization procedure for large |z| and small |z|, and similar resolutions 


were used in this case to generate the tabulated values. 


III. TESTING AND VERIFICATION OF RESULTS 


In addition to the previously published tables of spherical Bessel 
functions of complex argument previously mentioned [3], some tables of values 
for purely real and purely imaginary arguments were also available [11,12]. 
Excellent agreement was observed between values listed in these previously 
published tables and those obtained by the computer routines used to generate 


the tables of this report. 


In the case of purely imaginary arguments a connection can be made with 


modified spherical Bessel functions of real arguments using the relations 
[13]: 


i at®) * -\e1 Tat Yyf™) 
kx) “Vix Katy, 0) 


Ra -nmi/2 i/2 
Vx Bie If, Ge) se jy (xe" : ) 
Vz ae MY ee gett ty Gan) 
ntl of 
VK nt lf, () =F ac 1) ax ae Vp (x) - Dene Vp (x)]- (15) 


Reference 4 contains tables of the functions i,(x) and k(x) Equations 


(15) can readily be combined to give: 


nmi/2 


4, (4) =e 1 (x) 


-3(n+1)1i/2 


2 -n-1 
[1 ) - 5 (-1) 7k (x) - 


y, (ix) =e 


The values of j,(ix) and y, (ix) generated by using Eqs. (16) in 
bination with the tabulated values of 1,0) and k, (x) in Reference 4 were 
1 excellent agreement with the values obtained analyzing 5,¢2) and ¥,(2) with 


“wutines used to calculate the tables of this report. 


Although these concordances are encouraging, they still do not establish 
general validity of the results for complex arguments. As was discussed in 
Section II, the accuracy of the results obtained by using the recurrence 
relations decreases if the actual behavior of 4, (2) or Yn?) is contrary to 
the asymptotic behavior of these functions with respect to order. A method of 
testing the number of digits of accuracy of the resulting functions is 
available and consists simply of the substitution of the computed functions 
ifuto a cross product theorem. For the tables in this report, the following 


theorem was used for this purpose [14] 


§,(2) ¥,-1€2) ~ dy1(2) y,(e) = 2°. (7) 


The accuracy to which this theorem is satisfied guarantees that the 
values listed in the tables are accurate to six digits (with the sixth digit 


rounded ). 
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APPENDIX A 


Tables of 5,(2) and ¥n(2) for z between (1,1) and (10,10) 
and n between 0 and 20 


In the following tables (and those in Appendix B) the following notation 


is used: 


Two numbers enclosed in parentheses represent a complex 
number. For example, the expression (3,5) represents 
the complex number 3+5i. The notation Jn(z) represents 
the spherical Bessel function of the first kind (3,(2)] 
and Yn(z) represents the spherical Bessel function of 
the second kind fy, (2)]- Separate colum headings are 
provided for the real and imaginary parts of each of the 


functions. 
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APPENDIX B 


Tables of j,(z) and y,(z) for some extreme 


values of order and argument 


The following tables list values of 3,2) and ¥,(2) for certain extreme 
values of order and argument. Some question may arise as to the validity of 
the procedure outlined in Section II for computing functional values when the 
real and imaginary parts of the argument are of very aifferent order of 
magnitude (for example, when z = 103 + 10734). However, no difficulty arises 
even in such an extreme case. This is due to the fact that the recurrence 
relations were always initiated at a sufficiently large order m such that the 
real and imaginary parts of the largest desired order N of jn (2) were both 
seven orders of magnitude greater than the corresponding assumed starting 
values for in 2) This is all that is required for the iteration process to 
result in the required number of accurate digits. No limitations exist on the 
accuracy of the normalizing functions which were analyzed in such cases by 
expressions similar to Eq. (10). Once again, the infinite series 
representation of sinz and cosz were used for cases in which |z| is very 
small. Similarly, when the real part of z exceeded 21, the procedure 
discussed in connection with Eqs. (11) was used prior to invoking the 


appropriate FORTRAN functions. 
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